Installing and Loading Packages
- sf: import, manage and process vector-based geospatial data
- spatstat: perform point pattern analysis
- raster: read, write, manipulate, analyse and model of gridded spatial data
- maptools: provide tools for manipulating geographic data
- tmap: plot static point patterns maps or interactive maps
packages = c('sf', 'spatstat', 'raster', 'maptools', 'tmap')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}
Spatial Data Wrangling
Importing spatial data
st_read() of sf package is used to import the geospatial data
childcare_sf <- st_read("data/child-care-services-geojson.geojson") %>%
st_transform(crs=3414)
Reading layer `child-care-services-geojson' from data source
`C:\Users\user\Desktop\SMU\Y4S1\IS415\geniceee\IS415_blog\_posts\2021-09-15-hands-on-exercise-4\data\child-care-services-geojson.geojson'
using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84sg_sf <- st_read(dsn = "data", layer="CostalOutline")
Reading layer `CostalOutline' from data source
`C:\Users\user\Desktop\SMU\Y4S1\IS415\geniceee\IS415_blog\_posts\2021-09-15-hands-on-exercise-4\data'
using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21mpsz_sf <- st_read(dsn = "data",
layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source
`C:\Users\user\Desktop\SMU\Y4S1\IS415\geniceee\IS415_blog\_posts\2021-09-15-hands-on-exercise-4\data'
using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21Important to ensure that sg_sf and mpsz_sf are projected in the same projection system before using these data for analysis
print(st_crs(sg_sf))
Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["SVY21[WGS84]",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]print(st_crs(mpsz_sf))
Coordinate Reference System:
User input: SVY21
wkt:
PROJCRS["SVY21",
BASEGEOGCRS["SVY21[WGS84]",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]Although both sg_sf and mpsz_sf are projected in SVY21, the end of both prints indicate that the EPSG is 9001. This is a wrong EPSG code because the correct EPSG code for SVY21 should be 3414. There is therefore a need to assign the correct EPSG codes to sg_sf and mpsz_sf.
sg3414_sf <- st_set_crs(sg_sf, 3414)
mpsz3414_sf <- st_set_crs(mpsz_sf, 3414)
Mapping the geospatial datasets
It is useful to plot a map to show the geospatial data’s spatial patterns
tm_shape(sg_sf) +
tm_polygons() +
tm_shape(mpsz_sf) +
tm_polygons() +
tm_shape(childcare_sf)+
tm_dots()
Note: All the geospatial layers are within the same map extend. This shows that their referencing system and coordinate values are referring to similar spatial context
Alternatively, a pin map can be plotted. Its advantages are:
- Navigate and zoom around the map freely
- Query the information of each simple feature by clicking them
- Can change the background of the internet map layer. 3 are currently provided: ESRI.WorldGrayCanvas (default), OpenStreetMap and ESRI.WorldTopoMap
tmap_mode('view')
tm_shape(childcare_sf)+
tm_dots()
tmap_mode('plot')
Note: Always remember to switch back to plot mode afterwards because each interactive mode will consume a connection. Avoid displaying too many interactive maps in 1 RMarkdown document when publishing on Netlify.
Geospatial Data Wrangling
Many geospatial analysis packages require the input geospatial data in sp’s Spatial* classes
Converting sf data frames to sp’s Spatial* class
as_Spatial() of sf package is used to convert the data from sf data frame to sp’s Spatial* class
childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz3414_sf)
sg <- as_Spatial(sg3414_sf)
childcare
class : SpatialPointsDataFrame
features : 1545
extent : 11203.01, 45404.24, 25667.6, 49300.88 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables : 2
names : Name, Description
min values : kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>018989</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>1, MARINA BOULEVARD, #B1 - 01, ONE MARINA BOULEVARD, SINGAPORE 018989</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>THE LITTLE SKOOL-HOUSE INTERNATIONAL PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>08F73931F4A691F4</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
max values : kml_999, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>829646</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>200, PONGGOL SEVENTEENTH AVENUE, SINGAPORE 829646</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>RAFFLES KIDZ @ PUNGGOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>379D017BF244B0FA</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> mpsz
class : SpatialPolygonsDataFrame
features : 323
extent : 2667.538, 56396.44, 15748.72, 50256.33 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables : 15
names : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C, REGION_N, REGION_C, INC_CRC, FMEL_UPD_D, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area
min values : 1, 1, ADMIRALTY, AMSZ01, N, ANG MO KIO, AM, CENTRAL REGION, CR, 00F5E30B5C9B7AD8, 16409, 5092.8949, 19579.069, 871.554887798, 39437.9352703
max values : 323, 17, YUNNAN, YSSZ09, Y, YISHUN, YS, WEST REGION, WR, FFCCF172717C2EAF, 16409, 50424.7923, 49552.7904, 68083.9364708, 69748298.792 sg
class : SpatialPolygonsDataFrame
features : 60
extent : 2663.926, 56047.79, 16357.98, 50244.03 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables : 4
names : GDO_GID, MSLINK, MAPID, COSTAL_NAM
min values : 1, 1, 0, ISLAND LINK
max values : 60, 67, 0, SINGAPORE - MAIN ISLAND Converting the Spatial* class into generic sp format
spatstat requires the data to be in ppp object form. To convert Spatial* class to ppp object, the Spatial* class needs to be convert into Spatial object first.
childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")
childcare_sp
class : SpatialPoints
features : 1545
extent : 11203.01, 45404.24, 25667.6, 49300.88 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs sg_sp
class : SpatialPolygons
features : 60
extent : 2663.926, 56047.79, 16357.98, 50244.03 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs Converting generic sp format into spatstat’s ppp format
as.ppp()of spatstat is used to convert the Spatial object into spatstat’s ppp object format
childcare_ppp <- as(childcare_sp, "ppp")
childcare_ppp
Planar point pattern: 1545 points
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] unitsPlotting childcare_ppp
plot(childcare_ppp)
Summary statistics of childcare_ppp
summary(childcare_ppp)
Planar point pattern: 1545 points
Average intensity 1.91145e-06 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units
Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
(34200 x 23630 units)
Window area = 808287000 square unitsNote: In spatial point patterns analysis, the presence of duplicates is a significant issue. The statistical methodology used for spatial point patterns processes is based largely on the assumption that the process is simple (i.e. the points cannot be coincident)
Handling duplicated points
Check for any duplication in childcare_ppp
any(duplicated(childcare_ppp))
[1] TRUEmultiplicity() is used to count number of coincidence points. To know how many locations have >1 point event, refer to code chunk below:
sum(multiplicity(childcare_ppp) > 1)
[1] 128The output shows that there are 128 duplicated points. Plotting the data will allow us to view their locations
tmap_mode('view')
tm_shape(childcare) +
tm_dots(alpha=0.4,
size=0.05)
tmap_mode('plot')
To overcome this problem:
- Delete the duplicates (some useful point events will be lost)
- Jittering, which adds a small perturbation to the duplicate points so that they do not occupy the exact same space
- Make each point “unique” and then attach the duplicates of the points to the patterns as marks (as attributes of the points)
- Will however need analytical techniques that take into account these marks
Implementing Jittering Approach
childcare_ppp_jit <- rjitter(childcare_ppp,
retry=TRUE,
nsim=1,
drop=TRUE)
There are no duplicate points now
any(duplicated(childcare_ppp_jit))
[1] FALSECreating owin object
Good practice to confine the analysis within a geographical area when doing spatial point patterns analysis. In spatstat, an object called owin is specially designed to represent this polygonal region
sg_owin <- as(sg_sp, "owin")
plot(sg_owin)
Combining point events object and owin object
Extract childcare centres that are located within Singapore
childcareSG_ppp = childcare_ppp[sg_owin]
The output object combined both the point and polygon feature in 1 ppp object class
summary(childcareSG_ppp)
Planar point pattern: 1545 points
Average intensity 2.063463e-06 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units
Window: polygonal boundary
60 separate polygons (no holes)
vertices area relative.area
polygon 1 38 1.56140e+04 2.09e-05
polygon 2 735 4.69093e+06 6.27e-03
polygon 3 49 1.66986e+04 2.23e-05
polygon 4 76 3.12332e+05 4.17e-04
polygon 5 5141 6.36179e+08 8.50e-01
polygon 6 42 5.58317e+04 7.46e-05
polygon 7 67 1.31354e+06 1.75e-03
polygon 8 15 4.46420e+03 5.96e-06
polygon 9 14 5.46674e+03 7.30e-06
polygon 10 37 5.26194e+03 7.03e-06
polygon 11 53 3.44003e+04 4.59e-05
polygon 12 74 5.82234e+04 7.78e-05
polygon 13 69 5.63134e+04 7.52e-05
polygon 14 143 1.45139e+05 1.94e-04
polygon 15 165 3.38736e+05 4.52e-04
polygon 16 130 9.40465e+04 1.26e-04
polygon 17 19 1.80977e+03 2.42e-06
polygon 18 16 2.01046e+03 2.69e-06
polygon 19 93 4.30642e+05 5.75e-04
polygon 20 90 4.15092e+05 5.54e-04
polygon 21 721 1.92795e+06 2.57e-03
polygon 22 330 1.11896e+06 1.49e-03
polygon 23 115 9.28394e+05 1.24e-03
polygon 24 37 1.01705e+04 1.36e-05
polygon 25 25 1.66227e+04 2.22e-05
polygon 26 10 2.14507e+03 2.86e-06
polygon 27 190 2.02489e+05 2.70e-04
polygon 28 175 9.25904e+05 1.24e-03
polygon 29 1993 9.99217e+06 1.33e-02
polygon 30 38 2.42492e+04 3.24e-05
polygon 31 24 6.35239e+03 8.48e-06
polygon 32 53 6.35791e+05 8.49e-04
polygon 33 41 1.60161e+04 2.14e-05
polygon 34 22 2.54368e+03 3.40e-06
polygon 35 30 1.08382e+04 1.45e-05
polygon 36 327 2.16921e+06 2.90e-03
polygon 37 111 6.62927e+05 8.85e-04
polygon 38 90 1.15991e+05 1.55e-04
polygon 39 98 6.26829e+04 8.37e-05
polygon 40 415 3.25384e+06 4.35e-03
polygon 41 222 1.51142e+06 2.02e-03
polygon 42 107 6.33039e+05 8.45e-04
polygon 43 7 2.48299e+03 3.32e-06
polygon 44 17 3.28303e+04 4.38e-05
polygon 45 26 8.34758e+03 1.11e-05
polygon 46 177 4.67446e+05 6.24e-04
polygon 47 16 3.19460e+03 4.27e-06
polygon 48 15 4.87296e+03 6.51e-06
polygon 49 66 1.61841e+04 2.16e-05
polygon 50 149 5.63430e+06 7.53e-03
polygon 51 609 2.62570e+07 3.51e-02
polygon 52 8 7.82256e+03 1.04e-05
polygon 53 976 2.33447e+07 3.12e-02
polygon 54 55 8.25379e+04 1.10e-04
polygon 55 976 2.33447e+07 3.12e-02
polygon 56 61 3.33449e+05 4.45e-04
polygon 57 6 1.68410e+04 2.25e-05
polygon 58 4 9.45963e+03 1.26e-05
polygon 59 46 6.99702e+05 9.35e-04
polygon 60 13 7.00873e+04 9.36e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
(53380 x 33890 units)
Window area = 748741000 square units
Fraction of frame area: 0.414Plot the newly derived childcareSG_ppp
plot(childcareSG_ppp)
First-order Spatial Point Patterns Analysis
- Deriving KDE layer for visualising and exploring the intensity of point processes
- Performing Confirmatory Spatial Point Patterns Analysis by using Nearest Neighbour statistics
Kernel Density Estimation
Kernel density can be computed with density() of spatstat using the following configurations:
- bw.diggle() automatic bandwidth selection method. Other recommended methods are bw.CvL(), bw.scott() or bw.ppl()
- Smoothing kernel used is gaussian (default). Other smoothing methods are: “epanechnikov”, “quartic” or “disc”
- The intensity estimate is corrected for edge effect bias. The default is FALSE
kde_childcareSG_bw <- density(childcareSG_ppp,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
plot(kde_childcareSG_bw)
The output density values are way too small to comprehend. This is because the default unit of measurement of SVY21 is in m. As a result, the density values computed is in number of points/square metre
Retrieve bandwidth used to compute KDE layer
bw <- bw.diggle(childcareSG_ppp)
bw
sigma
298.4095 Rescaling KDE values
rescale() is used to covert the unit of measurement from m to km
childcareSG_ppp.km <- rescale(childcareSG_ppp, 1000, "km")
Re-run density() using the rescaled data set and plot the output KDE map
kde_childcareSG.bw <- density(childcareSG_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
plot(kde_childcareSG.bw)
Note: Output image looks identical to the earlier version, the only changes are in the data values
Working with different automatic bandwidth methods
Besides bw.diggle(), 3 other spatstat functions can be used to determine the bandwidth: bw.CvL(), bw.scott(), and bw.ppl().
Look at the bandwidth return by these automatic bandwidth calculation methods:
bw.CvL(childcareSG_ppp.km)
sigma
4.543278 bw.scott(childcareSG_ppp.km)
sigma.x sigma.y
2.224898 1.450966 bw.ppl(childcareSG_ppp.km)
sigma
0.3897114 bw.diggle(childcareSG_ppp.km)
sigma
0.2984095 bw.ppl() is suggested because it tends to produce more appropriate values when the pattern consists predominantly of tight clusters. However, if the purpose is to detect a single tight cluster in the midst of random noise, then bw.diggle() method seems to work best
kde_childcareSG.ppl <- density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")
6.1.3 Working with different kernel methods
By default, the kernel method used in density.ppp() is gaussian. There are three other options: Epanechnikov, Quartic and Dics
par(mfrow=c(2,2))
plot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="gaussian"),
main="Gaussian")
plot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="epanechnikov"),
main="Epanechnikov")
plot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="quartic"),
main="Quartic")
plot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="disc"),
main="Disc")
Fixed and Adaptive KDE
Computing KDE with fixed bandwidth
We will now compute a KDE layer by defining a bandwidth of 600m. Sigma value used is 0.6 because km is the unit of measurement for childcareSG_ppp.km
kde_childcareSG_600 <- density(childcareSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600)
Computing KDE with adpative bandwidth
Fixed bandwidth method is very sensitive to highly skew distribution of spatial point patterns over geographical units. Use adaptive bandwidth to overcome this issue
density.adaptive() of spatstat is used to derive adaptive kernel density estimation
kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)
Compare fixed and adaptive kernel density estimation outputs
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")
Converting KDE output into grid object
The result is the same, we just convert it so that it is suitable for mapping purposes
gridded_kde_childcareSG_bw <- as.SpatialGridDataFrame.im(kde_childcareSG.bw)
spplot(gridded_kde_childcareSG_bw)
Converting Gridded Output into Raster
raster() of raster package is used to convert the gridded kernel density objects into RasterLayer object
kde_childcareSG_bw_raster <- raster(gridded_kde_childcareSG_bw)
Properties of kde_childcareSG_bw_raster RasterLayer
Note that the CRS property is NA
kde_childcareSG_bw_raster
class : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348 (x, y)
extent : 2.663926, 56.04779, 16.35798, 50.24403 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : v
values : -8.476185e-15, 28.51831 (min, max)Assigning Projection Systems
CRS information will now be included in kde_childcareSG_bw_raster RasterLayer Note that the crs property is completed
projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster
class : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348 (x, y)
extent : 2.663926, 56.04779, 16.35798, 50.24403 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
source : memory
names : v
values : -8.476185e-15, 28.51831 (min, max)Visualising the output in tmap
tm_shape(kde_childcareSG_bw_raster) +
tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"), frame = FALSE)
Note: raster values are encoded explicitly onto the raster pixel using the values in “v” field
Comparing Spatial Point Patterns using KDE
KDE of childcare centres will be compared at Ponggol, Tampines, Chua Chu Kang and Jurong West planning areas
Extracting Study Area
pg = mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]
tm = mpsz[mpsz@data$PLN_AREA_N == "TAMPINES",]
ck = mpsz[mpsz@data$PLN_AREA_N == "CHOA CHU KANG",]
jw = mpsz[mpsz@data$PLN_AREA_N == "JURONG WEST",]
par(mfrow=c(2,2))
plot(pg, main = "Ponggol")
plot(tm, main = "Tampines")
plot(ck, main = "Choa Chu Kang")
plot(jw, main = "Jurong West")
Converting spatial point data frame into generic sp format
Convert SpatialPolygonsDataFrame layers into generic spatialpolygons layers
pg_sp = as(pg, "SpatialPolygons")
tm_sp = as(tm, "SpatialPolygons")
ck_sp = as(ck, "SpatialPolygons")
jw_sp = as(jw, "SpatialPolygons")
6.1.7.4 Creating Owin Object
Convert these SpatialPolygons objects into owin objects required by spatstat
pg_owin = as(pg_sp, "owin")
tm_owin = as(tm_sp, "owin")
ck_owin = as(ck_sp, "owin")
jw_owin = as(jw_sp, "owin")
Combining childcare points and the study area
Extract childcare centres within the regions for later analysis
childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]
rescale() is used to trasnform the unit of measurement from m to km
childcare_pg_ppp.km = rescale(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale(childcare_jw_ppp, 1000, "km")
par(mfrow=c(2,2))
plot(childcare_pg_ppp.km, main="Punggol")
plot(childcare_tm_ppp.km, main="Tampines")
plot(childcare_ck_ppp.km, main="Choa Chu Kang")
plot(childcare_jw_ppp.km, main="Jurong West")
Computing KDE
bw.diggle method is used to derive the bandwidth of each KDE
par(mfrow=c(2,2))
plot(density(childcare_pg_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Punggol")
plot(density(childcare_tm_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Tempines")
plot(density(childcare_ck_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Jurong West")
Computing fixed bandwidth KDE
par(mfrow=c(2,2))
plot(density(childcare_ck_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="Chou Chu Kang")
plot(density(childcare_jw_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="Jurong West")
plot(density(childcare_pg_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="Punggol")
plot(density(childcare_tm_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="Tampines")
Nearest Point Analysis
clarkevans.test() of statspat is used to perform Clark-Evans test of aggregation for a spatial point pattern
The test hypotheses are:
- Ho: distribution of childcare services are randomly distributed
- H1: distribution of childcare services are not randomly distributed
The 95% confident interval will be used
Testing spatial point patterns using Clark and Evans Test
clarkevans.test(childcareSG_ppp,
correction="none",
clipregion="sg_owin",
alternative=c("clustered"),
nsim=99)
Clark-Evans test
No edge correction
Monte Carlo test based on 99 simulations of CSR with fixed n
data: childcareSG_ppp
R = 0.54756, p-value = 0.01
alternative hypothesis: clustered (R < 1)Conclusion: Since p-value is less than 0.05, we can therefore reject the null hypothesis that distribution of childcare services are randomly distributed. Since the Nearest Neighbour index is < 1, we can infer that the distribution of childcare services resemble a clustered distribution
Clark and Evans Test: Choa Chu Kang planning area
clarkevans.test() of spatstat is used to performs Clark-Evans test of aggregation for childcare centres in Choa Chu Kang planning area
clarkevans.test(childcare_ck_ppp,
correction="none",
clipregion=NULL,
alternative=c("two.sided"),
nsim=999)
Clark-Evans test
No edge correction
Monte Carlo test based on 999 simulations of CSR with fixed n
data: childcare_ck_ppp
R = 0.93878, p-value = 0.09
alternative hypothesis: two-sidedClark and Evans Test: Tampines planning area
clarkevans.test(childcare_tm_ppp,
correction="none",
clipregion=NULL,
alternative=c("two.sided"),
nsim=999)
Clark-Evans test
No edge correction
Monte Carlo test based on 999 simulations of CSR with fixed n
data: childcare_tm_ppp
R = 0.76877, p-value = 0.002
alternative hypothesis: two-sidedSpatial Second-order Spatial Point Patterns Analysis
Analysing Spatial Point Process Using G-Function
G function measures the distribution of the distances from an arbitrary event to its nearest event. Gest() and envelope() of spatstat() package is used to compute G-function estimation and perform monte carlo simiulation respectively
Choa Chu Kang planning area
Computing G-Function Estimation
Gest() of spatstat package is used to compute G function
Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Choa Chu Kang are randomly distributed
- H1: distribution of childcare services at Choa Chu Kang are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
Monte Carlo test with G-function
G_CK.csr <- envelope(childcare_ck_ppp, Gest, nsim = 999)
Generating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60
.........70.........80.........90.........100.........110.........120
.........130.........140.........150.........160.........170.........180
.........190.........200.........210.........220.........230.........240
.........250.........260.........270.........280.........290.........300
.........310.........320.........330.........340.........350.........360
.........370.........380.........390.........400.........410.........420
.........430.........440.........450.........460.........470.........480
.........490.........500.........510.........520.........530.........540
.........550.........560.........570.........580.........590.........600
.........610.........620.........630.........640.........650.........660
.........670.........680.........690.........700.........710.........720
.........730.........740.........750.........760.........770.........780
.........790.........800.........810.........820.........830.........840
.........850.........860.........870.........880.........890.........900
.........910.........920.........930.........940.........950.........960
.........970.........980.........990........ 999.
Done.plot(G_CK.csr)
Tampines planning area
Computing G-Function Estimation
G_tm = Gest(childcare_tm_ppp, correction = "best")
plot(G_tm)
Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Tampines are randomly distributed
- H1: distribution of childcare services at Tampines are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
G_tm.csr <- envelope(childcare_tm_ppp, Gest, correction = "all", nsim = 999)
Generating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60
.........70.........80.........90.........100.........110.........120
.........130.........140.........150.........160.........170.........180
.........190.........200.........210.........220.........230.........240
.........250.........260.........270.........280.........290.........300
.........310.........320.........330.........340.........350.........360
.........370.........380.........390.........400.........410.........420
.........430.........440.........450.........460.........470.........480
.........490.........500.........510.........520.........530.........540
.........550.........560.........570.........580.........590.........600
.........610.........620.........630.........640.........650.........660
.........670.........680.........690.........700.........710.........720
.........730.........740.........750.........760.........770.........780
.........790.........800.........810.........820.........830.........840
.........850.........860.........870.........880.........890.........900
.........910.........920.........930.........940.........950.........960
.........970.........980.........990........ 999.
Done.plot(G_tm.csr)
Analysing Spatial Point Process Using F-Function
F function estimates the empty space function F(r) or its hazard rate h(r) from a point pattern in a window of arbitrary shape. Fest() and envelope() of spatstat() package is used to compute F-function estimation and perform monte carlo simiulation respectively
Choa Chu Kang planning area
Computing F-Function Estimation
Fest() of spatstat package is used to compute F-function
F_CK = Fest(childcare_ck_ppp)
plot(F_CK)
Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Choa Chu Kang are randomly distributed
- H1: distribution of childcare services at Choa Chu Kang are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
Monte Carlo test with F-function
F_CK.csr <- envelope(childcare_ck_ppp, Fest, nsim = 999)
Generating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60
.........70.........80.........90.........100.........110.........120
.........130.........140.........150.........160.........170.........180
.........190.........200.........210.........220.........230.........240
.........250.........260.........270.........280.........290.........300
.........310.........320.........330.........340.........350.........360
.........370.........380.........390.........400.........410.........420
.........430.........440.........450.........460.........470.........480
.........490.........500.........510.........520.........530.........540
.........550.........560.........570.........580.........590.........600
.........610.........620.........630.........640.........650.........660
.........670.........680.........690.........700.........710.........720
.........730.........740.........750.........760.........770.........780
.........790.........800.........810.........820.........830.........840
.........850.........860.........870.........880.........890.........900
.........910.........920.........930.........940.........950.........960
.........970.........980.........990........ 999.
Done.plot(F_CK.csr)
Tampines planning area
Computing F-Function Estimation
F_tm = Fest(childcare_tm_ppp, correction = "best")
plot(F_tm)
Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Tampines are randomly distributed
- H1: distribution of childcare services at Tampines are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
F_tm.csr <- envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)
Generating 999 simulations of CSR ...
1, 2, 3, ......10.........20.........30.........40.........50.........60
.........70.........80.........90.........100.........110.........120
.........130.........140.........150.........160.........170.........180
.........190.........200.........210.........220.........230.........240
.........250.........260.........270.........280.........290.........300
.........310.........320.........330.........340.........350.........360
.........370.........380.........390.........400.........410.........420
.........430.........440.........450.........460.........470.........480
.........490.........500.........510.........520.........530.........540
.........550.........560.........570.........580.........590.........600
.........610.........620.........630.........640.........650.........660
.........670.........680.........690.........700.........710.........720
.........730.........740.........750.........760.........770.........780
.........790.........800.........810.........820.........830.........840
.........850.........860.........870.........880.........890.........900
.........910.........920.........930.........940.........950.........960
.........970.........980.........990........ 999.
Done.plot(F_tm.csr)
Analysing Spatial Point Process Using K-Function
K-function measures the number of events found up to a given distance of any particular event. Kest() and envelope() of spatstat() package is used to compute K-function estimation and perform monte carlo simiulation respectively
Choa Chu Kang planning area
Computing K-Function Estimate
K_ck = Kest(childcare_ck_ppp, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(m)")
Performing Complete Random Spatial Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Choa Chu Kang are randomly distributed
- H1: distribution of childcare services at Choa Chu Kang are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
K_ck.csr <- envelope(childcare_ck_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.plot(K_ck.csr, . - r ~ r, xlab="d", ylab="K(d)-r")
Tampines planning area
Computing K-Function Estimation
K_tm = Kest(childcare_tm_ppp, correction = "Ripley")
plot(K_tm, . -r ~ r,
ylab= "K(d)-r", xlab = "d(m)",
xlim=c(0,1000))
Performing Complete Random Spatial Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Tampines are randomly distributed
- H1: distribution of childcare services at Tampines are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
K_tm.csr <- envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.
Analysing Spatial Point Process Using L-Function
Lest() and envelope() of spatstat package is used to compute L-function estimation and perform monte carlo simulation test respectively
Choa Chu Kang planning area
####Computing L-Function Estimation
L_ck = Lest(childcare_ck_ppp, correction = "Ripley")
plot(L_ck, . -r ~ r,
ylab= "L(d)-r", xlab = "d(m)")
Performing Complete Spatial Random Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Choa Chu Kang are randomly distributed
- H1: The distribution of childcare services at Choa Chu Kang are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
L_ck.csr <- envelope(childcare_ck_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")
Tampines planning area
Computing L-Function Estimate
L_tm = Lest(childcare_tm_ppp, correction = "Ripley")
plot(L_tm, . -r ~ r,
ylab= "L(d)-r", xlab = "d(m)",
xlim=c(0,1000))
Performing Complete Spatial Random Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted:
- Ho: distribution of childcare services at Tampines are randomly distributed
- H1: distribution of childcare services at Tampines are not randomly distributed
The null hypothesis will be rejected if p-value < 0.001 (alpha value)
L_tm.csr <- envelope(childcare_tm_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.